convert.to.m.i<-function(i,s.list,player,u.tbl) {
  s<-s.list[i]
  u.tbl.sub<-subset(u.tbl,u.tbl[,player]==s)
  
  u.tbl.sub[,(player+2)]
}

convert.to.m<-function(u.tbl) {
  s1.list<-levels(as.factor(u.tbl[,1]))
  s2.list<-levels(as.factor(u.tbl[,2]))
  
  pi1<-t(sapply(1:length(s1.list),convert.to.m.i,s.list=s1.list,1,u.tbl))
  pi2<-t(sapply(1:length(s2.list),convert.to.m.i,s.list=s2.list,2,u.tbl))
  
  list(pi1,pi2)
}

quantal.response<-function(gamma,s.list,s.dist,pi.m) {
  prob.m<-matrix(rep(s.dist[,2],length(s.list)),nrow=length(s.list),byrow=TRUE)
  sc<-exp(gamma*(apply(pi.m*prob.m,1,sum)))
  cbind(s.list,sc/sum(sc))
}

find.qre<-function(gamma,u.tbl,tol=1e-4,max.iter=100,mode="tbl") {
  gamma1<-gamma[1]
  gamma2<-gamma[2]
  
  if(mode=="tbl") {
    pi.list<-convert.to.m(u.tbl)
    s1.list<-as.numeric(levels(as.factor(u.tbl[,1])))
    s2.list<-as.numeric(levels(as.factor(u.tbl[,2])))
  }
  if(mode=="m") {
    pi.list<-u.tbl
    s1.list<-as.numeric(dimnames(pi.list[[1]])[[1]])
    s2.list<-as.numeric(dimnames(pi.list[[2]])[[1]])
  }
  
#  browser()
  pi1<-pi.list[[1]]
  pi2<-pi.list[[2]]
  
  sc1<-exp(gamma1*(apply(pi1,1,mean)))    #each row a strategy
  sc2<-exp(gamma2*(apply(pi2,1,mean)))
  qr1<-cbind(s1.list,sc1/sum(sc1))
  qr2<-cbind(s2.list,sc2/sum(sc2))
  
  acc<-10000
  counter<-0
  
  qr.list<-NULL
  
  while((acc>tol)&(counter<max.iter)) {
    qr1.new<-quantal.response(gamma1,s1.list,qr2,pi1)
    qr2.new<-quantal.response(gamma2,s2.list,qr1,pi2)
    
    #    acc<-sum(abs(qr1.new-qr1))+sum(abs(qr2.new-qr2))
    
    acc<-sum(qr1[,2]*(log(qr1[,2])-log(qr1.new[,2])))+sum(qr2[,2]*(log(qr2[,2])-log(qr2.new[,2])))
    
#    if(is.na(acc)) browser()
    
    qr1<-qr1.new
    qr2<-qr2.new
    
    counter<-counter+1
    qr.list<-rbind(qr.list,c(counter,acc))
  }
  
  # calculate payoff
  
  pm1<-matrix(qr1[,2],nrow=length(s1.list))%*%matrix(qr2[,2],ncol=length(s2.list))
  pm2<-matrix(qr2[,2],nrow=length(s2.list))%*%matrix(qr1[,2],ncol=length(s1.list))
  
#  browser()
  
  list(qr1,qr2,qr.list,sum(pm1*pi1),sum(pm2*pi2))  
}

